If one doesn't change any parameters or initial conditions will initialization still run on every call to solve? (I'm trying to understand what kind of overhead this might add when solving SDEs / jump problems where one needs to run many identical simulations to calculate statistics.)

For clarity, the semantics that we're settling towards is the following:

1. Always default to running initialization on all problem types. This is because you can have parameter dependencies and that kind of stuff in any situation. This system is somewhat required for the generality MTK allows, otherwise you can end up in some odd incorrect cases. And with that, at solve we validate that we are solving the problem we think we're solving.
2. In many cases this initialization is trivial. If it is trivial, then the nonlinear solve has u0 === nothing. The initialization system just does a check for if initializationdata.u0 === nothing and skips any solving in that case. Note that trivial initialization is not just cases for which there are no initialization equations but also the case where MTK is able to analytically handle the equations (for example, linear equations with tearing).
3. If we are in the trivial case, then we still may need to resolve the explicit system which is the u0map, i.e. the map from the initialization parameters to that of the full system. This is done at the solve time. Because verification takes longer than running such a function, we verify by simply running the explicit map. This ensures correctness. You can bypass that with NoInit()... but you're really splitting hairs there because it's an O(n) function called one time so it's pretty much guaranteed to be less than 5% of the cost of one ODE solver step or nonlinear solver step... so if the solver does more than one step you've already made that trivial.
4. For simplicity, we will then also evaluate that map during remake when trivial. So basically, remake will get a keyword argument lazy = initializationdata.u0 === nothing, so you can explicitly turn that on/off at will. If you turn it on for a non-trivial case then you'd have to incur the nonlinear solve cost (and you have to using NonlinearSolve or some nonlinear solver, with an error message for help), and that's a requirement that cannot be relaxed for obvious reasons. But at least the trivial case does not require any solver and it's super cheap so we can run it without any issue. The only downside is that it might give someone the wrong idea that it's always okay to assume that the values from the remake will be non-lazy.

That makes the easy cases easy, the hard cases possible, and with essentially zero measurable overhead in any case I've tried.

Requires:

• feat: add infrastructure for initialization of different problem types SciMLBase.jl#885
• feat: add optional history function to ProblemState SymbolicIndexingInterface.jl#109
• feat: support initialize_dae! for SDEIntegrator StochasticDiffEq.jl#594
• feat: add initialization support DelayDiffEq.jl#301
• feat: add initialization support StochasticDelayDiffEq.jl#84

DAEProblem is not included in the added tests because it needs some special treatment (#3240). I need to find ways to re-create the tested scenarios for JumpProblem. NonlinearProblem/NonlinearLeastSquaresProblem will also be tested once I get the infrastructure in NonlinearSolve.

Requires SciML/NonlinearSolve.jl#517 for NonlinearSolve

Working on adding tests/support for initialization of SCCNonlinearProblem. There's also some fairly significant performance bugs I've found in this infrastructure which I'll address.

MTK tests pass locally. Waiting on DelayDiffEq.jl and NonlinearSolve.jl (via SciML/NonlinearSolve.jl#518) to be tagged to push the changes. The performance issue mentioned above has also been fixed locally.

This is good to go. InterfaceI failures are due to OrdinaryDiffEq, InterfaceII failures are fixed in #3284.